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ABSTRACT 

There are two distinct regimes of the first order Fermi acceleration at shocks. The first is a hnear (test 
particle) regime in which most of the shock energy goes into thermal and bulk motion of the plasma. 
The second is an efficient regime when it goes into accelerated particles. Although the transition region 
between them is narrow, we identify the factors that drive the system to a self- organized critical state 
between those two. Using an analytic solution, we determine this critical state and calculate the spectra 
and maximum energy of accelerated particles. 

Subject headings: acceleration of particles — cosmic rays — diffusion — hydrodynamics — radiation 
mechanisms mon-thcrmal — shock waves 



1. INTRODUCTION 

It is now generally recognized that most of the observed 
gamma radiation derives in one way or another from ac- 
celerated particles. Radio and x-ray spectral components 
from a variety of astrophysical objects are believed to have 
a similar origin. High energy neutrinos, whose detection 
is on the program for the future and existing water/ice 



dctootora, muat alao bo related to the ultra high energy 
cosmic Y&ya (LTHECR) . Their origin is, in turn, a mystorjr 



and the huge Auger detector complex is now being built 
to elucidate it ( [Blandford 1999] , [Cronin 1999| ). 

There has been essential progress in our understand- 
ing of how accelerated particles produce the radiation de- 
tected. The models concentrate on the synchrotron and 
inverse Compton emission for electrons and on the 7 ray 
and neutrino production in pp and P7 reactions for pro- 
tons. However, the primary spectrum of accelerated parti- 
cles remains a stumbling block making predictions of oth- 
erwise similar models so differe nt (e.g., Mannheim et al., 
1999| , [Waxman fc BahcaU 1999| ). 

The "standard" mechanism of particle acceleration, ca- 
pable of producing nonthermal power-law spectra extend- 
ing over many decades in energy is the I-order Fermi or 
diffusive shock acceleration. It was originally suggested 
to explain the origin of galactic cosmic rays (CRs). For 
the purposes of the high-energy radiation and UHECR it 
is usually adopted as an axiom, and mostly only in its 
simplest, test particle (TP) or linear realization. In par- 
ticular, it is assumed that any strong nonrelativistic shock 
routinely produces a E~'^ spectrum of protons and/or elec- 
trons. In fact this spectrum arises from a simple formula 
F ^ £'-(''+2)/(»'-i) where r is the shock compression (for 
r = 4). However, it is valid only if the shock thickness 
is much smaller than the particle mean free path. This, 
in turn, is true only if the energy content of accelerated 
particles is small compared to the shock energy (inefficient 
acceleration) so that the shock structure is maintained by 
the thermal, not by the high energy particles. Otherwise, 



the accelerated particles create the shock structure on their 
own and if so, then obviously on a scale that is larger or of 
the order of their mean free path, thus making the above 
formula invalid. Therefore, the TP regime requires a very 
low CR number density (the rate of injection v into the 
acceleration process), which appears to be impossible in 
the parameter ran ge of inter est. It has been inf erred from 
observations (e.g., Lee 1982), s imulations (e .g., Bennet & 
EUison 199"5| ), and theory (e.g., |A-Ialkov 199^ ) that the CR 
number density ncR at a strong shock must be ~ 10~^ 
of the background density ni upstream. It is important 
to emphasize here that when the actual injection rate v 
exceeds the critical value (denote it 1^2), the test parti- 
cle {E~'^) solution simply does not exist. Simple mea- 
sures, such as calculating corrections, are intrinsically in- 
adequate. What happens is that two other solutions with 
considerably higher efficiencies branch off at a somewhat 
lower injection rate i/i < ^"2, one of which disappears again 
at v — V2^ together with the test particle solution. 

Thus, it seems to be difficult to put an accelerating 
shock into a regime in which the CR energy production 
rate (acceleration efficiency) could be gradually adjusted 
by changing parameters. It is either too low (TP regime) 
or it is close to unity. Note, that this situation is quite 
suggestive of that occurring in phase transitions or bifur- 
cations. 

Generally, neither of those extreme regimes provide an 
adequate description of particle spectra and related emis- 
sion. Nevertheless, we argue in this Letter that despite this 
apparent lack of regulation ability, shocks must be still ca- 
pable of self-regulation and self- organization. The transi- 
tion region between the two acceleration regimes (critical 
region) is very narrow in control parameters like u. On the 
other hand, the self-regulation can work efficiently only 
when the parameters are within this region. This require- 
ment determines them, and also resolves the question of 
the mechanism of self-regulation. The above consideration 
is similar to the concept of self- organized criticality (SOC) 



(e.g., iBak Tang fc Wiesenfeld 1987| , |Hwa fc Kardar 1992 , 
Diamond fc Hahm 19951)! 



2. FORMULATION OF THE PROBLEM 



We use the diffusion-convection equation (e.g., [Drury 



1983 ) for describing the distribution of high energy par- 
We assume that the gaseous discontinu- 

and 
Thus, 



tides (CRs 

ity (also called the subshock) is located at x 

the shock propagates in the positive x- direction 

the flow velocity in the shock frame can be represented 

as V{x) = —u{x) where the (positive) flow speed u{x) 

jumps from U2 = u(0— ) downstream to uq = u(O-l-) > U2 

across the subshock and then gradually increases up to 

ui = u(-l-oo) > uo- In a steady state the equation reads 






IdM dj_ 
3 dx dp' 



(1) 



where f{x,p) is the isotropic (in the local fluid frame) part 
of the particle distribution. This is assumed to vanish far 
upstream (/ — > 0, a; — > 00), while the only bounded so- 
lution downstream is obviously f{x,p) = foip) = /(0,p). 
The most plausible assumption about the cosmic ray dif- 
fusivity k{p) is that of the Bohm type, i.e., k(p) = 
Kp^ / yjl + p'^ (the particle momentum p is normalized to 



mc). In other words k scales as the gyroradius, k 



;{?)■ 



The reference diffusivity K depends on the SB/B level 
of the MHD turbulence that scatters the particles in pitch 



angle. The minimum value for K would be K 



VeS if 



6B ^ B. Note that this plain parameterization of this im- 
portant quantity is perhaps the most serious incomplete- 
ness of the theory which will be discussed later. 

To include the backreaction of accelerated particles on 
the plasma flow three further equations are needed. First, 
one simply writes the conservation of the momentum flux 
in the smooth part of the shock transition (a; > 0, the 
so-called CR-precursor) 



Pr 



pu 



Piul, 



x>0 



where Pc is the pressure of the high energy particles 



Prix) — — mc / 

3 7^,, 



p'^dp 

Vp^ + i 



f{p,x) 



(2) 



(3) 



It is assumed here that there are no particles with mo- 
menta p > pi (they leave the shock vicinity because there 
are no MHD waves with sufficiently long wave length A, 
since the cyclotron resonance requires p ~ A). The mo- 
mentum region < p < po cannot be described by equa- 
tion {m and the behavior of f{p) a,t p ^ p o is described 
by the injection parameters po and /(po) (Malkov 1997, 
[M97]). The plasma density p{x) can be eliminated from 
equation (0) by using the continuity equation pu — piui. 
Finally, the subshock strength Vg can be expressed throug h 



the Mach number M at a; = 00 (e.g.. Landau & Lifshitz) 



"2 



7 + 1 



7 - 1 + 2i?T+iAf- 



(4) 



where the precursor compression R = ui/uq and 7 is the 
adiabatic index of the plasma. 

The system of equations (0,00) describes in a self- 
consistent manner the particle spectrum and the flow 



structure. An efficient way to solve it is to reduce this 
system to one integral equation (M97). A key dependent 
variable is an integral transform of the flow profile u{x) 
with a kernel suggested by an asymptotic solution of the 
system (||)-(||) which has the form 



f{x,p) ^ /o(p)exp 



6K 



where 



*= / u{x')dx' 
/o 



and the spectral index downstream q{p) 
The integral transform is as follows 



-dlnfo/dlnp. 



U{p) = — I cxp 
wi Jq- 



3k{p) 
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and it is related to q{p) through the following formula 



lip) = 



3 



dlnU 



dlnp rsRU{p) 



(6) 



Thus, once U{p) is found both the flow proflle and the par- 
ticle distribution can be determined by inverting transform 
dq) and integrating equation (o) . Now, using the linearity 
of equation (0) {pu = const) , we derive the integral equa- 
tion for U by applying the transformation (pf) to the x— 
derivative of equation (0) (M97). The result reads 



U(t) 
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K{t') K{t)q{t) 
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U{t") 



(7) 



where t = Inp, to.i = Inpo.i- Here the injection parameter 
47r mc^ 



- - ^5 2P0MP0) 

3 piuf 

is related to R by means of the following equation 
1^ = Kpo (1 - R-^) 



(8) 
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The equations (||0ia) form a closed system that can be 
easily solved numerically. We analyze the results in the 
next section. 

3. MECHANISMS OF CRITICAL SELF-ORGANIZATION 

The critical nature of this acceleration process is best 
seen in variables R, v. The quantity i? — 1 is a measure of 
shock modiflcation produced by CRs, in fact (i? — \)/R = 
Pc{Q)/ piu\ (eq.||]) and may be regarded as an order pa- 
rameter. The injection rate v characterizes the CR density 
at the shock front and can be tentatively treated as a con- 
trol parameter. It is convenient to plot the function v{R) 
instead of R{i') (using equation 0), since R{v) is not al- 
ways a single-valued function. Fig. |l|. 



The injection rate v at the subshock should be calcu- 
lated given Tsiji) (M97) with the self-consistent determi- 
nation of the flow compression R on the basis of the Riv) 
dependence obtained. However, in view of its critical char- 
acter, this solution can be physically meaningful only in 
regimes far from criticality, i.e., when Rk, \ (test particle 
regime) or i? ::^ 1 (efficient acceleration). But, it is diffi- 
cult to see how this system could stably evolve remaining 
in one of these two regimes. Indeed, if v is subcritical it 
will inevitably become supercritical when pi is sufficiently 
high. Once it happened, however, the strong subshock re- 
duction (equation Q) will reduce v and drive the system 
back to the critical regime. Fig. 0. 

The maximum momentum pi is subject to self- 
regulation as well. Indeed, when R^ \, the generation 
and propagation of Alfven waves is characterized by strong 
inclination of the characteristics of wave transport equa- 
tion towards larger wavenumbers k on the k — x plane 
due to wave compression. Thus, considering particles with 
p <, Pi inside the precursor, one sees that they are in reso- 
nance with waves that must have been excited by particles 
with p > pi further upstream but, there are no particles 
with p > pi. Therefore, the required waves can be excited 
only locally by the same particles with p <, pi which sub- 
stantially diminishes the amplitude of waves that are in 
resonance with particles from the interval Pi/R < p < pi- 
(The left inequality arises from the resonance condition 
kcp ~ eB/mc and the frequency conservation along the 
characteristics ku{x) = const). This will worsen the con- 
finement of these particles to the shock front. The quan- 
titative study of this process is the subject of current re- 
search. What can be inferred from Fig. g now, is that 
the decrease of pi straightens and rise the curve v{R), so 
that it returns to the monotonic behaviour. However once 
the actual injection becomes subcritical (and thus R^ 1) 
then pi will grow again restoring the two extrema on the 
curve u{R). 

The above dilemma is quite typical for dynamical sys- 
tems that are close to criticality or marginal stability. A 
natural way to resolve it consists in collapsing the extrema 
into an inflection point so that a self- organized critical 
(SOC) acceleration regime is established being determined 
by the conditions iy'{R*) — v"{R*) = 0. These condi- 
tions not only determine unique critical values R* and 
V* = i'{R*) but also yield the maximum momentum pi 
as a function of M, which is shown in Fig. g. A few 
particle spectra that develop in the SOC states for differ- 
ent Mach numbers are shown in Fig. H, along with the 
asymptotic M = oo non SOC spectrum. The latter can 
be calculated in a closed form (M97). Note that the hard- 
ening of the spectra in about the last decade below the 
cut-off is entirely due to the abrupt cut-off itself. 

4. DISCUSSION 

The detailed microphysics behind the SOC is extremely 
complex and must include the self-consistent turbulence 
evolution and particle acceleration with their strong back- 
reaction on the shock structure. We have simplified it and 
argued that the most important dynamical components, 
the bulk plasma flow and the high-energy particles, must 
be in a balance that constitutes a certain equipartition 
of the shock energy between the two. This was done by 
identifying the factors that prevent either of them from 



prevailing alone. 

The above situation is similar to that in e.g., a simple 
sandpile paradigm of the SOC. It is impossible (in fact un- 
necessary) to describe the individual grain dynamics, but 
it is clear that when the critical macroscopic characteristic 
of the system (the slope of the sandpile) becomes too steep 
due to the action of external factors, like tilting of the en- 
tire system or addition of sand at the top, the sandpile 
relaxes bringing the slope to its critical magnitude. 

5. POSSIBLE FEEDBACK FROM OBSERVATIONS 

Perhaps the most significant observational aspect is the 
particle spectrum. Although the conversion of detected 
radiation spectra into the primary particle spectra is am- 
biguous, in some cases it may be compared with the theory. 
The most striking prediction is that in shocks with very 
high maximum particle energy, the spectra must be harder 
than g = 4 (or 2 in the normalization f{E)dE). This is 
because of a very low injection requirement for efficient 
acceleration in such shocks (Fig.|l]) . If we (conservatively) 
set ncn/ni ~ 10~^, then v ^ lO^^cpo/mul which may 
easily exceed 1^2 (local maximum on the J^(i?) curve) al- 
ready for pi ^ 10^ — 10^, putting the acceleration into 
a strongly nonlinear regime. This should have important 
observational consequences. 

First, the particle energy is concentrated at the upper 
cut-off instead of being evenly distributed over the loga- 
rithmic energy bands as in the test particle E~'^ solution. 
This makes the upper bounds on CR gen erat ed neutrino 
fiuxes (see e.g., Bahcall & Waxman 1999 and Mannheim 
et at, 199£ and references therein) rather ambiguous. In- 



deed, the UHECR spectrum is normalized to the observed 
one at i?norm ~ 10^^ eV while being obscured by the galac- 
tic background at the energies E <^ 10^^ eV. Thus, CR 
spectra harder than E~^ imply tha t the upper limit on 
the n eutrino fluxes e.g., derived by Waxman fc Bahcalj 
1999| should be even lower than the E~"' spectrum implies 
for the energies E^y < 5- 10^"^ eV {E,,/Ecr ^ 0.03 is a typi- 
cal energy relation). According to the same logic, it should 
be increased for higher energies. Note that the upper cut- 
offs in individual shocks contributing to the UHECR must 
be still much higher than -Enorm to validate our simpli- 
fied handling of particle losses. On the other hand if the 
sources with i^inax < ^-norm Contribute significantly, the 
measurements at -Bnorm tell us nothing about their nor- 
malizations and the upper bound on neutrino fluxes may 
be increased up to the level dictated by the C R observa- 
tions in the lower energy range, as suggested by Mannheinj 
et al., 1999| . This scenario is supported by Fig. ||, provided 
that there are many strong shocks in the ensemble. The 
observed CR steep power-law spectrum is then essentially 
a superposition of flatter or even non-power-law spectra 
from individual sources properly distributed in E'max- 

To summarize our conclusions, the main factor that 
should determine the particle primary spectra, and thus 
the neutrino flux, is how the accelerating shocks are dis- 
tributed in cut-off momenta, which in a SOC state means 
in Mach numbers. There is no universal spectral form for 
individual shocks at the current state of the theory (except 
a not quite representative case of M — > 00). Therefore, 
one should understand particle losses mechanisms, since 
they determine the shock structure and thus the spectra. 



directly and through E^i^^. These mechanisms are insepa- 
rable from the dynamics of strong compressible MHD tur- 
bulence generated by those same particles. The further 
progress in its study will improve our understanding of 
the acceleration process and related radiation. 
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Fig. 1. — The nonlinear response of an accelerating shock (characterized by the precursor compression R) to the thermal injection f 
given in the form of the function ulR) calculated for the fixed injection momentum po = 10~^, Mach number M = 150 and for different 
pi = 100; 550; 10*; 10^; lO". The critical value (see text) p* = 550. The TP regime is limited to the region R~l. 

Fig. 2. — Bifurcation diagram corresponding to the set of response curves shown in fig. hi. Since u and pi are in reality dynamic rather than 
control parameters the response curve moves towards the bifurcation curve drawn with the heavy line. It corresponds to the curve marked 
by Pi = 550 in Fig 
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Fig. 3. — Maximum momentum pi versus Mach number M calculated in the SOC states (see text). 

Fig. 4. — The spectral indices in the SOC states shown for different Mach numbers and the corresponding maximum momenta from Fig. 
bI (solid curves). For comparison, the dashed curve shows the asymptotic case M = oo, pi = 10^^ and u = 4: ■ 10~^. For larger, e.g., SOC 
injection values this spectrum would have been cut in mildly relativistic region. 



